/*************************************************************************
Purpose: Standard libraries are loaded, and commonly used types (such as
         structures) and functions are defined here.
**************************************************************************
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
#ifndef GAIA_DEFINITIONS
#define GAIA_DEFINITIONS
#define GAIA_VERSION 0.620

//Un-comment next line to compile with GSL. Comment if GSL ISN'T installed.
//Remember to change LFLAGS in Makefile so it doesn't link ANY GSL libs.
//#define GSL_HERE

//Load C libraries which are common to nearly all GAIA programs:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>//For realpath().
//Allows the use of cos() and sin(), among other functions:
#include <math.h>
//Used for searching directories:
#include <dirent.h>
//Used for sizing notification intervals in loops on the fly and
//for providing estimates of time remaining:
#include <time.h>
//The following are used for fork().
#include <errno.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/wait.h>

//Load C++ libraries which are common to nearly all GAIA programs:

#include <iostream>
//The following libraries allow setw() to set widths in cout statements:
#include <iomanip>
#include <fstream>
#include <string>
#include <vector>
//Used for the random rand() function:
#include <cstdlib>
//Used for the min() function:
#include <algorithm>
//Used for string streams:
#include <sstream>

//Load custom library of LOTS of mathematical constants.
#include "double_constants.hpp"

//GSL libraries.
#ifdef GSL_HERE
  #include <gsl/gsl_linalg.h>
  #include <gsl/gsl_blas.h>
  #include <gsl/gsl_sf_erf.h>
  #include <gsl/gsl_sf_legendre.h>
  #include <gsl/gsl_errno.h>
  #include <gsl/gsl_fft_real.h>
  #include <gsl/gsl_fft_halfcomplex.h>
  #include <gsl/gsl_rng.h>
  #include <gsl/gsl_randist.h>
  #include <gsl/gsl_cdf.h>
#endif

//Used to load FES2004 routines- note that fes.h is not actually in this
//folder as "fes.h" would seem to indicate.  It's actually in the library
//that FES's Makefile generates and is included in the g++ compile statement
//in the Makefile in this folder.
#include "fes.h"

//So that I don't have to precede vector and string with std::
using namespace std;

//////////////////////////////////////////////////////////////////////////
//List of definitions and physical constants.
//////////////////////////////////////////////////////////////////////////

#define ACCEL_READ_FORMAT "%lld %lf %lf %le %d"
#define ACCEL_WRITE_FORMAT "%11lld%10.4f%10.4f%13.5e%12d\n"
#define ACCEL_HEADER "   UTC Time  Latitude Longitude Accel(m/s^2)\n"
#define ORBIT_READ_FORMAT   "%lld %lf %lf %le %d %le %le %le %le %le %le"
#define ORBIT_WRITE_FORMAT  "%11lld%10.4f%10.4f%13.5e%12d%21.7e%21.7e%21.7e%21.7e%21.7e%21.7e\n"
#define ORBIT_HEADER "   UTC Time  Latitude Longitude Accel(m/s^2)   Ascending                   xA                   xB                   yA                   yB                   zA                   zB\n"
#define MAP_HEADER_UNSTRUCTURED "     Latitude     Longitude              Output   Elevation(m)\n"
#define DOUBLE_FORMAT "%20.9f\n"
#define CUSTOM_GRID_LS_SUBFOLDER "custom-grid"

//My constants for speeding up various unit conversions.
static const double seconds_in_year = 3600*24*365.2422;//Average days/yr.
static const double inverse_seconds_in_year = 1.0/seconds_in_year;
static const double rad2deg = 180.0/pi;
static const double deg2rad = pi/180.0;

//Any added periods also need to go here: is_this_period_recognized() and
//convert_title_to_parameter().

//Tidal periods, in seconds.  Taken from FES source code file prediction.c
//and output from my program fes_tidal_periods.cpp.
static const double q1_period =          26.8683567132*3600.0;
static const double o1_period =          25.8193417018*3600.0;
static const double k1_period =          23.9344695923*3600.0;
static const double two_n2_period =      12.9053744727*3600.0;
static const double n2_period =          12.6583482359*3600.0;
static const double m2_period =          12.4206012119*3600.0;
static const double k2_period =          11.9672347962*3600.0;
static const double s2_period =          12.0000000000*3600.0;
static const double p1_period =          24.0658902255*3600.0;
static const double m4_period =           6.2103006102*3600.0;
static const double s1_period =          24.0000000000*3600.0;
static const double mf_period =         327.8589689145*3600.0;
static const double mm_period =         661.3092048547*3600.0;
static const double mtm_period =        219.1903862122*3600.0;
static const double msqm_period =       170.2990185857*3600.0;

/* Replaced by FES values on 2009-05-15
//Tidal periods, in seconds.  Taken from Shailen's thesis, except
//for annual, semi-annual which are based on 365.2422 days/year.
static const double q1_period =                26.868356705*3600.0;
static const double o1_period =                25.819341660*3600.0;
static const double k1_period =                23.934469223*3600.0;
static const double two_n2_period =            12.904800000*3600.0;//Wunsch
static const double n2_period =                12.658348083*3600.0;
static const double m2_period =                12.420600891*3600.0;
static const double k2_period =                11.967234612*3600.0;
static const double s2_period =                12.000000000*3600.0;
static const double p1_period =                24.065889359*3600.0;
static const double m4_period =                m2_period/2.0;
static const double mf_period =               327.859008789*3600.0;

//From tide glossary (downloaded).
static const double mm_period =               661.309204855*3600.0;

static const double mtm_period =              219.120000000*3600.0;//Wunsch
static const double msqm_period =             170.304000000*3600.0;//Wunsch
//*/

//m2 f + 1cyc/yr
static const double m2_plus_period =           12.403026588*3600.0;
//m2 f - 1cyc/yr
static const double m2_minus_period =          12.438225067*3600.0;

static const double lunar_sidereal_period =   655.728027344*3600.0;

static const double semiannual_period =       seconds_in_year/2.0;
static const double annual_period =           seconds_in_year;
static const double four_year_period =        seconds_in_year*4;

//From Ray and Luthcke 2006
static const double k1_alias_period =          2725.4*24*3600.0;
static const double k2_alias_period =          1362.7*24*3600.0;
static const double o1_alias_period =          13.6*24*3600.0;
static const double p1_alias_period =          171.2*24*3600.0;

//Tidal angular frequencies, in radians/second.
static const double m2_omega = twoPi/(m2_period);
static const double m4_omega = twoPi/(m4_period);
static const double m2_plus_omega = twoPi/(m2_plus_period);
static const double m2_minus_omega = twoPi/(m2_minus_period);
static const double k1_omega = twoPi/(k1_period);
static const double s1_omega = twoPi/(s1_period);
static const double s2_omega = twoPi/(s2_period);
static const double o1_omega = twoPi/(o1_period);
static const double p1_omega = twoPi/(p1_period);
static const double n2_omega = twoPi/(n2_period);
static const double two_n2_omega = twoPi/(two_n2_period);
static const double k2_omega = twoPi/(k2_period);
static const double msqm_omega = twoPi/(msqm_period);
static const double mtm_omega = twoPi/(mtm_period);
static const double mf_omega = twoPi/(mf_period);
static const double mm_omega = twoPi/(mm_period);
static const double q1_omega = twoPi/(q1_period);
static const double lunar_sidereal_omega = twoPi/(lunar_sidereal_period);
static const double semiannual_omega = twoPi/(semiannual_period);
static const double annual_omega = twoPi/(annual_period);
static const double four_year_omega = twoPi/(four_year_period);

static const double k1_alias_omega = twoPi/(k1_alias_period);
static const double k2_alias_omega = twoPi/(k2_alias_period);
static const double p1_alias_omega = twoPi/(p1_alias_period);
static const double o1_alias_omega = twoPi/(o1_alias_period);

//Whenever a time value is stored in UTC seconds past J2000.0 but is used
//to calculate FES, this value needs to be subtracted from it so that the
//times are relative to FES's epoch time of 1950-01-01 00:00:00.0.
//fes_utc_ext_epoch is used for times that are sent to fes_access().
static const long long fes_utc_ext_epoch = -1577880000LL;//1950-1-1.0 rel to J2000
//This FES epoch value was reasoned as follows:
//1.  FES values seem to be in Universal Time relative to 01-JAN-1950 00:00:00.000 according to the FES SDK (I call this UT1950).
//2.  char2sec.e (used to obtain mult_times from tide gauge files) converts character date:time values to time in seconds past J2000.
//3.  Thus, I ran the following conversion:
//bryan@Prometheus:~/GAIA/software/time_conversion$ ./char2sec.e 01-JAN-1950 00:00:00.000
//and got the answer: -1577880000.0000000
//This value is "UT1950 - J2000".
//UPDATE: On 2009-08-19, John Wahr noticed an epoch of 1900-01-01 in prediction.c:
//fes_utc_int_epoch is used for phases that are relative to internal 1900 epoch.
static const long long fes_utc_int_epoch = -3155716800LL;//1900-1-1.0 rel to J2000

//Converts Julian centuries to Julian days.
static const double jcy2jday = 36525.0;

//Converts seconds to Julian centuries.
static const double sec2jcy = 1.0/(86400.0*36525.0);

//Converts seconds to Julian days.
static const double sec2jday = 1.0/86400.0;

//If two numbers are less than "TOLERANCE" apart, they are "identical".
static const double tolerance = 1.0E-8;

//Some applications require a higher tolerance.
static const double high_tolerance = 1.0E-3;

//Many loading functions require use c-style strings, which need a maximum
//number of characters to be defined.
//2009-09-09 Upgrade: New grid.wcover/mean_depth needs 400 from 300.
static const int max_length = 400;

//FES loading code requires a c-style array buffer with length FES_BUFFER
#define FES_BUFFER 30

//Topography code complains if grid pt is farther from topo pt than:
static const double max_topo_res = 50.0;

//Quadratic mean radius of the earth (in km) from Wikipedia.
static const double earth_avg_radius = 6372.795477598;
//Polar radius of the earth (in km) from gd2xyz.e.
static const double earth_polar_radius = 6356.7513016;
//Average earth gravity acceleration, m/s^2.
static const double earth_little_g = 9.80665;

//Programs use planet_x_radius instead so other planets can be analyzed.
static const double planet_avg_radius   = earth_avg_radius;
static const double planet_polar_radius = earth_polar_radius;
static const double little_g = earth_little_g;

//static const double GEOID_r_major=6378.1370;//WGS84 semi-major axis (km)
static const double GEOID_r_major=6378.136;//from ./gd2xyz.e -H (km)
//static const double GEOID_r_minor=6356.752314245;//WGS84 semi-minor axis (km)
static const double GEOID_r_minor=6356.751301555;//derived from ./gd2xyz.e flattening (km)
static const double GEOID_f=1.0-GEOID_r_minor/GEOID_r_major;//WGS84 flattening
static const double GEOID_e2=2*GEOID_f-GEOID_f*GEOID_f;//WGS84 eccentricity squared
static const double GEOID_ed2=GEOID_r_major*GEOID_r_major/(GEOID_r_minor*GEOID_r_minor)-1.0;//WGS84 eccentricity derived
static const double GEOID_e=sqrt(GEOID_e2);//WGS84 eccentricity

//Degrees to km arclength.
static const double deg2km = deg2rad*planet_avg_radius;
//km arclength to degrees.
static const double km2deg = rad2deg/planet_avg_radius;

//Newton's gravitational constant
static const double Big_G = 6.67428E-11;

//For right now, uniform, non-time-dependent density (kg/m^3).
static const double sea_density = 1035.0;

//Convert from psia pressure to cm of water.
//P in Pascals = density*g*(cm/100)
//cm = Pascals*100/(density*g)
//1.0 psi = 6894.75729 Pascals (google)
//cm = 100*6894.75729*psi/(sea_density*little_g)
static const double psi2cm = 100*6894.75729/(sea_density*little_g);

//////////////////////////////////////////////////////////////////////////
//List of custom data types.
//////////////////////////////////////////////////////////////////////////

typedef struct {

  string basefolder;//The folder (relative to GAIA/) that has input file.
  string filename;//input filename.
  string output_filename;//Occasionally, input data is saved using this.
  string output_folder;//Occasionally, input data is saved using this.
  int noload;//Set to 1 if the input isn't supposed to be loaded.
  int nocache;//Set to 1 if input.binindices isn't to be loaded from cache.
  int no_orbit_data;//Set to 1 if the input file has no orbit data.
  string cachefolder;//Folder where cached inputs are stored.
  long long segmentsize;//When loading input in segments, size of each segment.

  //Sometimes only some input data is used for inversion- other pts 4 verify
  int split_input;
  long long chunk,buffer;//Used in split_input.
} input_options_s;

typedef struct {

  //Each run of this program can produce multiple plots, so the data
  //is stored in a vector of vectors and the titles and units are also
  //stored in vectors.
  vector<string> titles;//Title that ends up on IDL plots.
  vector<string> x_units,y_units;
  vector< vector<string> > legends;//One for each plot's lines.
  //Vector of vectors of doubles for the output XY points.
  //1st index is the plot number- each plot is in a different file and pic.
  //2nd index is the dataseries- each plot can have multiple lines.
  //3rd index is the index of the x/y value in each dataseries.
  vector< vector< vector<double> > > x_values,y_values;

  //Used to interpolate arbitrary xy timeseries to input time series.
  vector<long long> closest_indices;
  vector<double> differences;//for closest_indices, arbitrary - input.
  double numindices;
  double max_diff, min_diff;
  long long smallest_index, largest_index;
  vector<int> parameter_type;//copied from results.
  vector<double> parameter;//copied from results.
  long long current_parameter;
  string cachefolder;//Folder where preloaded arbitrary time series are stored.
} xy_s;

typedef struct {

  vector<double> lat,lon,accel; //Coords of midpoint and relative accel.
  vector<double> xA,xB,yA,yB,zA,zB;//Coords of GRACE A,B in xyz system.

  vector<long long> time;//GPS time.
  //In order to convert these time values to the same coordinate system as
  //the FES model (UT 1950):
  long long fes_epoch;
  //In order to convert these time values to the same coordinate system as
  //any new model- Richard Ray's, Padman's, etc.:
  long long newmodel_epoch;

  vector<short> ascending;//1(0) if satellites are headed north(south).
  vector<long long> binindices;//Used w region code- depends on input and grid
  //Translates binindices into the indices of each region's grid, not the
  //indices in the global grid as with binindices.
  vector<long long> regional_binindices;

  //When data is loaded for a region, store the indices of each input time
  //(i.e. where that input time/accel/position can be found in the main
  //input file) in ipir.  This is mainly useful for up continuation, but
  //originally it was used for debugging region code.
  vector<long long> indices_of_points_in_region;//debug of region code!
  vector<long long> indices_of_points_in_region2;//debug of region code!

  long long first_time;//First time in input file.
  long long last_time;//Last time in input file.
  long long epoch;//Currently set to the midpoint.  Subtracted from all times.

  //This vector is the same length as binindices and is either 0 (pt used in
  //inversion) or 1 (pt not used in inversion, only for verification).
  vector<long long> verification_points;

  //Controls compare_input_files().
  int file_comparison_method;

  //Used for computing upward continuation of just ocean tides, sep. nontides, or land tides.
  vector<double> tides_accel,nontides_accel,landtides_accel;
  
  vector<long long> bad_indices;//In combine_(), pts with big km_error are deleted.

  //Holds interpolated versions of each arbitrary parameter's time series.
  xy_s arbitrary;

  //Holds copy of options so I don't have to pass input_options around.
  input_options_s options;
} input_s;

typedef struct {

  int type;//Controls how topography data is loaded.
  string basefolder;//Folder where topography data is stored.
  string filename;//Filename of topography data.
  vector<double> lat,lon,elev;

  //Used for SRTM30+ which is equal-grid and stored as 2-byte ints.
  vector< vector <short> > hr_elev;
  //hr_lats/lons are much shorter than hr_elev; they only store each row/column.
  vector<double> hr_lats, hr_lons;

  //1 - List of lats,lons,elevs.
  //2 - Elevs in 2D equal-grid grid.
  int format;
} topo_s;

typedef struct {

  string basefolder; //Folder where temp data is stored.
  string filename; //Filename of temp data.
  vector<double> lat,lon,temp;
} temp_s;

typedef struct {

  topo_s topo;
  temp_s temp;
} ancillary_s;

typedef struct {

  int type;
  double lat_spacing, lon_spacing, cap_size;
  double boundaries[4];
  string cachefolder;//Folder where cached grids are stored.
  int noload;//Set to 1 if you don't want the grid loaded from cache.
  int area_type;//Specifies area of disks, squares, etc.

  int support_type;//Controls how supporting grid is formed.
  int support_lat_multiplier,support_lon_multiplier;

  //Needed because automatic cap sizing alters these.
  double original_lat_spacing, original_lon_spacing, original_cap_size;

  //Grid type 11, region type 3, analysis_list_choice 1 all use this point.
  double specified_lat,specified_lon;

  //In kilometers. When snapping grid/region to specified pt, the snap
  //is considered successful if the mismatch is less than this distance.
  double specified_tolerance;

  //Useful for cases when region code wraps around grid code.
  string base;//Usually either "grid" or "region".

  int topo_type;//Copy of ancillary.topo.type for grid cache substring.

  int latlon;//0 disables latlon, >1 for equal-grid latlon grid (see create_options in gaia.cpp).
  int current_latlon;//Is used to tell define_grid to make a latlon grid, NOT latlon!
  long long latlon_lat_multiplier, latlon_lon_multiplier;// >= 1, makes latlon denser than grid.
  
  //Used only in region code.
  int parallel;//Switches on/off region parallel processing code.
  double overlap_distance;//Used to control how many "edge" points there are.

  int analysis_list_cache_only;//Skips any regions not in ls cache.
  int analysis_list_choice;//Used to determine which regions are analysed.
  //manual_list holds region indices that are used in some a_list_choices.
  vector<long long> manual_list;
} grid_options_s;

typedef struct {

  string filename;
  double mask;//Mask values occur on land.
  //Used for loading FES amplitudes and phases on 2D grid- obsolete/redundant?
  vector< vector<double> > amplitudes, phases;
  //Used for generic multi-parameter output on 2D grid.
  vector< vector< vector<double> > > outputs;
  vector<double> lat,lon;//Defines coords of 2D grid.
  vector< vector<double> > elev;//Elevation (m), probably that of nearby unstructured grid point.
  int discard_redundant_values;//Zero keeps them, anything else gets rid of redundant (lon=360) values in FES.
} latlon_s;

typedef struct {

  //These variables are used by all grid types.
  vector<double> lat, lon, elev, x, y, z, distance_to_land, distance_to_ocean, wcover, wdepth;
  //Keeps a short, efficient list of the different lat bands in the grid.
  vector<double> different_lats;
  vector<long long> chosen_for_region;//debug!
  vector<double> area;//Area of each grid point for kg->cm conversion
  vector<double> temp;//Time-avg temp of each gr pt for kg->cm conversion(C)
  vector<double> density;//Density of each grid point for kg->cm conversion
  vector<double> cm2kg;//Used for cm-> kg conversion.
  vector<double> radius;//Radius(km) of each gr pt for disk mascons
  vector<double> radius_angle;//Radius(radians) for disk mascons

  //These variables are used to calculate the boundaries of bins
  //when the analysis method chosen uses bins (this is nearly always the case).
  double maxlat,minlat,maxlon,minlon;//Boundaries.
  double minlonsmall,maxlonsmall,minlonlarge,maxlonlarge;//Alt. boundaries.
  vector<double> lat_halfwidth;//Defined in terms of degrees or km.
  vector<double> lon_halfwidth;//Defined in terms of degrees or km.

  //Holds coordinates for supporting grid points.  1st index is main gr pt
  //that is "supported" by all the grid points in the 2nd index.
  //support_weights holds weighting factors for support grid points.
  vector<vector<double> > support_lat, support_lon, support_x, support_y, support_z, support_weights;

  //main_weights holds weighting factors for main grid points.
  vector<double> main_weights;

  //Every time the grid is rotated (in any direction) this variable
  //(in degrees) is used to perform the rotation. It's a temp var only.
  double rotate;

  //Holds copy of options so I don't have to pass grid_options around.
  grid_options_s options;
  
  //Holds equal-grid latlon grids.
  latlon_s latlon;

  //Only used by region code.

  //Contains all input data in this region.
  input_s input;

  //Each pt in regional_grid has an index to the original grid structure.
  //Some pts in regional_grid are "center" pts and some are "edge" pts.
  //First index is region #, second indices run over the grid points
  //in that region.
  vector< vector<long long> > all_grid_indices, center_grid_indices, edge_grid_indices;

  //Created by build_analysis_region_list().
  //Controls which regions are analyzed in main().
  vector<long long> indices_to_analyze;
  
  //Used to match FES pts with mascons.
  vector< vector<long long> > fes_lat_indices, fes_lon_indices;
  vector< vector< vector<long long> > > fes_index;
  vector<double> fes_areas;//Area for FES pts at different lats.
} grid_s;

typedef struct {

  int type;//Determines which model to load.
  string basefolder; //The folder (relative to GAIA/) that has model files.
  vector<string> filenames; //Model filenames.
  string cachefolder;//Folder where cached models are stored.
  int noload;//Set to 1 if you don't want the model loaded from cache.

  double halfwidth;//Controls interpolations using Gaussian smoothing.
  double averaging_radius;//Controls averaging.

  //If 1, load_model saves values within avg_rad to results in a
  //similar manner to type 1003.
  int record_nearby_values;

  //Controls how tide gauge time series are constructed.
  int tg_type;
} model_options_s;

typedef struct {

  string title;//Title of this model, "FES 2004", "R. Ray", "Padman", etc.
  grid_s grid; //Allows model to be loaded onto its own unique grid.

  //1st index is M2, K1, etc.  2nd index is model.grid index.
  vector< vector<double> > amplitudes, phases, sin, cos;
  //Stores "M2", "K1", etc to label first index of "amplitudes", "phases".
  vector<string> titles;

  //The 1st index in amplitudes, phases refers to tide period.  This tide
  //period is stored here, in seconds.
  vector<double> periods;

  //The 1st index in amplitudes, phases refers to tide period.  The angular
  //frequency ("omega") corresponding to this period is stored here, in
  //radians per second.
  vector<double> omegas;

  string basefolder;//Folder (relative to GAIA/) that has model files.
  vector<string> filenames; //Model filenames.
  string cachefolder;//Folder where cached models are stored.

  //When comparing the model to tide gauge stations, this vector of vectors
  //holds the model's predicted sealevel height at each station.
  //1st index is the tide gauge station.
  //2nd index is the time index.
  vector< vector<double> > sealevels;//(Mean has been removed!)

  //Used if record_nearby_values==1.
  double avg_amp,avg_phase,closest_amp,closest_phase,closest_dist;

  //Holds copy of options so I don't have to pass model_options around.
  model_options_s options;
} model_s;

typedef struct {

  //Allows args and phase to be calculated using:
  //1 - Shailen's thesis
  //2 - FES code
  //3 - Wenzel/Simon 1994
  int type;

  double input_period;//Period in question.
  //That input_period is then recognized and associated with:
  vector<int> numbers;//6 Doodson numbers, stored from 1-6 (0 is empty!).

  double tau;//Mean lunar time (UNITS?!??).
  double s;//Mean lunar lon (deg).
  double h;//Mean solar lon (deg).
  double p;//Mean lunar perigee lon (deg).
  double np;//Negative mean lon of ascending lunar node (deg).
  double ps;//Mean solar perigee lon (deg).
  vector<double> args;//All of the above args for easy looping.
  vector<long double> omegas;//In rad/sec for the above args, #6 IS VERY SMALL!
  vector<double> periods;//In seconds for the above omegas.
  double output_period;//Calced from Doodson #s/omegas, should = input_period.

  //Used in eq 4.1,2 of Shailen's thesis, values from Table 2.2 on p25.
  int h2kj_sign;

  //6 Doodson titles, stored from 1-6 (0 is empty!).
  vector<string> short_titles, long_titles;

  double phase;//Eq. 2.31 from Shailen's thesis, omega*t + beta
  double time;//Used to calculate phase; in seconds rel. to epoch.
} doodson_s;

typedef struct {

  int recognized;//Was current period recognized as a period defined here?
  int tide;//1 if this period is a tide (ex: M2) or not (ex: annual).
  int parameter_type;//Used in convert_title_to_parameter- for const, linear.
  string title; //Title of this period, if recognized.
  double period;//Period in seconds.
  double omega;//twoPi/(period in seconds) - plugs into cos(omega*t).
  int output_type;//Just like any entry in results.output_type.
  int sincos;//If period is encoded as sincos, this is 1. amp/phase = 0

  doodson_s doodson;//If tide, contains AT LEAST Doodson nums. Args/phase?
} period_info_s;

typedef struct {

  string basefolder; //Folder where tide gauge data is stored.

  //The following vector holds the titles of each tide gauge station.
  vector<string> titles;//Tide gauge desc "Abbas's Aasiaat data"

  //Some tide gauge data only contains amplitudes and phases.  Because
  //of this fact, each tide gauge station needs to have have a "1" in this
  //vector (meaning it already has filled its entries in sealevels)
  //or a "0", meaning that it only has amplitudes and phases.
  vector<int> includes_timeseries;//1 if station includes time series else 0
  //AP data stored here.
  vector< vector<double> > periods, phases, amplitudes;
  vector< vector<int> > tide_booleans;
  vector< vector<string> > period_titles;

  //Grid points are simply coords, elevs of the stations.
  grid_s grid;
  grid_options_s grid_options;

  //Vector of vectors of long longs and doubles for the time series data.
  //1st index is the tide gauge station.
  //2nd index is the time index.
  vector< vector<long long> >   times;//UTC seconds past J2000.0
  vector< vector<double> > sealevels;//(Mean has been removed!)
  vector< vector<double> > airpressures;//in mbars or hectoPascals.
  vector< vector<double> > salinities;

  //Once inversion is interpolated to in situ stations, use amps/phases to
  //predict the RESIDUAL sealevel at all the times for which tide gauge data
  //exists.(Mean has been removed!)
  vector< vector<double> > input_sealevels;

  //When using native FES 2004 routines to predict sealevel height at these
  //stations, it's often useful to store the Julian dates that FES expects-
  //these are relative to 1950-01-01 00:00:00.
  vector< vector<double> > fes_julian_times;

  //The model is evaluated NEAR each station, and results are stored here.
  //This copy of model has model.grid.lat.size() == # of stations.
  model_s model;

  //This copy of model_options exists solely so that its type can be
  //set independently of the user-set model_options.type.  This way, even if
  //the user chooses to load the model's amps/phases, the sea level height
  //at each station, at each time can be loaded inside load_ancillary().
  model_options_s model_options;

  //Each tide gauge has a list of grid indices which are "nearby"
  //They're used to weight the mascon residual signal to get a guess
  //at what the mascon residual signal is right at the tg station.
  vector< vector<long long> > nearby_grid_indices;
  vector< vector<double> > nearby_grid_distances;
  vector< vector<double> > weighting_factors;

  //The resulting guess is stored in this vector of vectors, in either
  //amplitude/phase or cosamp/sinamp form, depending on the value in
  //results.output_type.
  vector< vector<double> > output_values;
} tg_s;

typedef struct {

  tg_s tg;
  int comparison_type;//Controls how model is compared to in situ data.
  string cachefolder;//Folder where preloaded in situ data is stored.
  int noload;//Set to 1 if you don't want in situ data loaded from cache.
  //1 to load, 0 to skip.
  int loadUHSLC,loadPSMSL,loadABBAS,loadAP,loadCEC,loadCanada,loadDART,loadFRAM;

  int init_only;//1/0 tells load_insitu to only load locations/load everything.

  string name1,name2;
  vector< vector<double> > series1,series2;

  //Used to create synth time series.
  long long synth_interval;//Between data points, in seconds.
  long long num_synth_points;
  long long synth_offset;//In seconds, usually relative to J2000. Defines time of first data point.
} insitu_s;

typedef struct {
  //Requires GSL libraries.
  #ifdef GSL_HERE
    //This is the first type of "H" matrix storage I programmed, on 2008-01-13.
    gsl_matrix * gsl;
  #endif
  //Unless I make an assignment statement to h_matrix that's not in the
  //switch-case, I get a warning on Azonic.  Stupid!
  int stop_uninitialized_warning;
} h_matrix_s;

typedef struct {
  //Requires GSL libraries.
  #ifdef GSL_HERE
    gsl_matrix * A;//H^T*H
    gsl_vector * b;//H^T*y
    long long dimension;//Dimension of A (square, symmetric matrix) and b.

    gsl_matrix * R;//Regularization matrix- same size as A.
    gsl_matrix * F;//Intermediate step to making R- same size as A.
  #endif
  string reg_title;//Contains name of the regularization method used in R.

  //Same as GSL above but faster population.
  vector <double> stl_b,stl_x;
  vector < vector<double> > stl_A;
} least_squares_s;

typedef struct {

  //The following variables are used for partial save/load code.

  //This variable is 0 if no partial_save_details are found in a subfolder
  //which is being loaded (indicating that the folder is completed).
  //It's 1 if partial_save_details files are found, or if a new folder
  //is being created (b/c it's always a good idea to make partial saves).
  int partial_save_code_used;

  //This variable holds the number of regions that are completed already
  //AND have been requested by the current region.indices_to_analyze.
  //(Second one is similar but for partially completed requested regions.)
  long long completed_requested_regions,partial_requested_regions;

  //All vectors are init'd to have an entry for each requested region, in
  //requested order.

  //These are filenames of requested partial save details files.
  vector <string> filenames;

  //0 = no detail file found, no region cache file found. Start from scratch
  //1 = no detail file found, reg cache file found. Presumed completed by old code.
  //2 = detail file found, not finished. Load region's cache file and start
  //    populating from the first unused input point.
  //3 = detail file found, finished. Just load from cache and finish!
  //    (The same thing should happen if an old ls cache is loaded, one from
  //    before details files. But this should only happen if NO details
  //    files were found at all- otherwise this region was never calc'ed!)
  vector <int> initial_completeness;

  vector<long long> input_numpoints;

  //This was the last index summed over in the partial save cache file.
  vector<long long> previous_input_index;

  //Yes, these should all be the same number.
  vector<long long> input_epoch;

  //For each region; epoch has already been added back.
  vector<long long> first_input_time,last_input_time;

  //If a region is complete, same as ps_first_input_time,ps_last_input_time.
  vector<long long> first_used_input_time,last_used_input_time;
} ps_s;

typedef struct {

  int type;
  //Value of parameter (unimportant for constant, secular trend,
  //but holds the period in seconds for sine/cosine terms).
  vector<double> parameter;

  //Int codes that indicate what kind of parameter is being passed.
  //See create_options() in gaia.cpp for more info.
  vector<int> parameter_type;

  //Describes format of output, controls write_output(), etc.
  // 1 - Unstructured map.
  // 2 - 2D plots- each plot can have multiple curves and a legend.
  // 3 - Map separate points (used in support grid maps).
  // 4 - Maps of separate points with labels.
  // 5 - Map using 2D latlon equal-grid format.
  int output_choice;

  //0 for all parameter types other than sinusoids.  0 for all non-tidal
  //periods (ex: annual).  1 for periods that are tides. (ex: M2).
  vector<int> parameter_tide;

  //Analysis types that work in mascon mass amplitude space require this:
  int h_matrix_choice;//Controls how create_h_matrix creates the H matrix.

  //Disables loading the least squares matrices and vectors A,b from cache.
  //NOTE: Does NOT disable saving to cache!
  int least_squares_noload;

  //Folder where cache subfolders reside.  Loaded from options.txt
  string least_squares_cachefolder;

  //Folder where cache FILES reside, calced in ls_cache_init.  Relative.
  string least_squares_cachesubfolder;

  //List of files in subfolder with filenames that have more than 2 chars.
  vector<string> least_squares_cache_filenames;

  //List of region indices that ls_cache_filenames refers to.
  vector<long long> least_squares_cache_region_indices;

  //Because create_least_squares_in_parallel is run multiple times in a parallel loop,
  //it's a good idea to decide whether or not the function will save TO
  //the cache in subfolder least_squares_cachesubfolder or load FROM that
  //folder.  The subfolder name also needs to be decided before the loop!
  //If least_squares_cache_load == 0, fctn should not affect cache at all.
  //Technically, this is true.  But it's not possible for ls_cache_init()
  //to return a value of 0- it pretty much always sets it to -1,1,2!
  //This is probably a good thing- why would I not want to save to cache?
  //If least_squares_cache_load == 1, fctn should LOAD from subfolder using
  //"exact parameter match" simple loading code.
  //If least_squares_cache_load == 2, fctn should LOAD from subfolder using
  //"reordered or subset" complex loading code.
  //If least_squares_cache_load == -1, fctn should SAVE to subfolder.
  //NOTE: The only effect least_squares_noload has is to force any values
  //of 1,2 to -1.  It only has this effect in
  //least_squares_cache_initialize(), so this forcing doesn't take place
  //when the value would be 0 anyway.
  int least_squares_cache_load;

  //After I wrote ls_cache_load, I realized there's a different type of
  //subset matching I might want to do - grid_subset matching (regions
  //that have a smaller overlap distance than cached region's matrices).
  //This variable is set to 1 if this is happening, and it also
  //forces ls_cache_load to 2.  I believe that even exact parameter
  //matches can be handled by the reordered parameter loading code, so this
  //should be okay.
  int grid_subset;

  //If least_squares_cache_load == 2, need to use this vector to load
  //the A,b matrices from the files in least_squares_cachesubfolder.  The
  //first (0'th) index in this vector tells the loading code where to find
  //the 0'th currently desired parameter in the cached parameter list.
  vector<long long> least_squares_parameter_indices;

  //If/when a least squares cache folder is found, record the value used
  //to create its A matrix and b vector- h_size (may be different later!)
  long long ls_cache_h_size;
  //Also need parameter_type:
  vector<int> ls_cache_parameter_type;

  //Copy of the exact region struct formed (rather, loaded, hopefully) by
  //region_options in the least squares cache.  This can be slightly
  //different from desired region struct (because I'm allowing loading
  //of least squares subfolders that were made with HIGHER overlap
  //distances.
  grid_s ls_region;

  //Copy of region currently desired, in its entirety.  Center and edge
  //grid indices will be compared to those in ls_region for grid subsets.
  //current_grid holds a copy of the ACTUAL COMPLETE grid used in main(),
  //which must be exactly the same for cache grid and desired grid.
  //region.
  grid_s current_region, current_grid;

  //Sometimes, it's useful to allow perform_analysis() to know which region
  //is currently being used.  The current region index is stored here.
  long long region_index;

  //Sometimes, it's useful to let least_squares routines know which REQUESTED region
  //is currently being used.  The current REQUESTED region index is stored here.
  //In other words, if there are 6 regions but region.indices_to_analyze = [5,4,1]
  //then the first loop would have region_index = 5 and requested_index = 0,
  //the second loop would have region_index = 4 and requested_index = 1,
  //and the third loop would have region_index = 1 and requested_index = 2.
  //requested_index is most often used to access ps.* vectors because they're the
  //same size as region.indices_to_analyze.
  int requested_index;

  //Analysis types that return sin and cos coefficients have this variable
  //set to "1", types that return amp and phase have a setting of "0".
  int sincos;

  //Usually set in perform_analysis() to control regularize() results.
  //0 - Gaussian smoothing matrix.
  //1 - Tikhonov
  int regularization_choice;

  //Parameters used in regularization()- ex: gaussian radius
  double regparam1,regparam2;
  int regparamint1;

  //Controls whether mascons are verified by upward continuing them.
  //0 - don't verify.
  //1 - verify.
  int upward_continuation_verification_choice;
  
  //If != 0, any verification choice calling for upward continuation will
  //skip that process and instead load the subfolder "folder_for_skipping".
  int skip_upward_continuation_during_verification;
  string folder_for_skipping;

  int numprocesses;//How many processes run in parallel (region or time)

  //Controls whether create_least_squares_in_parallel... runs in parallel.
  //0 - Do not run in parallel (default).
  //1 - Run in parallel using MPI - use same number of processes as region.
  int time_parallel;

  //Controls the vector over which the grid loop in create_least_squares()
  //cycles through grid points.  Thus, grid points "too far" away from
  //input pt don't have to be calculated in the contribution to the A,b
  //matrices.
  //0 - Cycle over all gr pts by creating trivial zone_of_influence v o vs.
  //1 - Spherically symmetric zone of radius zone_radius.
  int zone_of_influence_choice;

  double zone_radius;

  vector< vector<long long> > zone_of_influence;

  //Trivia about zone_of_influence size:
  double zone_avg_pts;
  long long zone_max_pts,zone_min_pts;

  //For disk mascons, Legendre polynomials are expanded from n=0 to n=nmax.
  long long nmax;
  int uses_nmax;//Non zero only if hmc uses nmax, so ls cache init knows.

  //For disk mascons, each disk has this radius:
  double disk_radius;
  int uses_disk_radius;//!=0 only if hmc uses radius, so ls cache init knows.

  int uses_support_grid;//So ls cache init knows to look at supp gr params

  //Regularization requires incrementing the A matrix with lambda*R, where
  //R is a regularization matrix from regularize().
  double lambda;

  //Controls refined inversion loop.
  double convergence;

  //Controls time series with input filters.
  double cutoff_freq;

  //Controls plotting FFTs where independent variable is period.
  double min_period,max_period;

  string prefix;//When using region code, this is set to "Region   1", etc.

  int found_file;//Used in least_squares_load_from_cache().

  //This structure contains variables used for partial save/loading code.
  ps_s ps;

  //HMC 202 populates using STL multivec, moves data to ls cache.
  int already_moved_to_cache;

  //Even if STL is called for, each region may load completed cache files,
  //in which case GSL containers are required. BUT RESET FOR EACH REGION!
  //Set to 0 at beginning of clsip(), set to 1 if COMPLETED region found.
  int loading_completed_matrices;

  //Documents current phase format of results.
  //0 - Unknown.
  //1 - Relative to J2000.
  //2 - Greenwich phase lags.
  //3 - Relative to input.epoch.
  int current_phase_format;

  //Controls adjust_phases() which alters current_phase_format.
  //See adjust_phases() in functions.cpp for more info.
  int adjust_phases_choice;

  //If 0, copy_regional_results only copies center grid pts to master gr.
  //If 1, copy_regional_results copies all grid points to master grid.
  int copy_overlap_points;

  //1 - Writes analysis results (also in overlap grid points!) to disk.
  //This variable is set by the user in create_options. The other variable
  //called copy_overlap_points is set to 0 all the time except in
  //overlap_results where it is set to 1.
  int save_overlap_points;

  //1 - Saves all the input data found in each region.
  int save_regional_input_data;

  //1 - Immediately prints results at first grid point.
  int print_immediate_results;

  //Like a_o.parameter_type but 2 entries for sinusoids.  Values:
  //0    - Constant (cm)
  //1    - Secular trend (cm/year)
  //2    - Acceleration (cm/year^2)
  //3-99 - Time^n term
  //101  - Cosine coefficient (cm)
  //102  - Sine coefficient (cm)
  //103  - Amplitude (cm)
  //104  - Phase (degrees)
  vector<int> output_type;

  //Like a_o.parameter_tide but has 2 entries for sinusoids.  Values:
  //0 for all parameter types other than sinusoids.
  //0 for all non-tidal periods (ex: annual).
  //1 for periods that are tides. (ex: M2).
  vector<int> output_tide;

  //1 if this analysis type uses least squares routines, 0 otherwise.
  int uses_least_squares;

  //Allows perform_analysis() to be run only to redefine param vectors.
  int inside_load_options;
} analysis_options_s;

typedef struct {

  //The angular frequency ("omega") corresponding to every period is
  //stored here, in radians per second.
  vector<double> omegas;

  //Each run of this program can produce multiple plots, so the data
  //is stored in a vector of vectors and the titles and units are also
  //stored in vectors.
  vector<string> titles;//Title that ends up on IDL plots.
  vector<string> units;//Units that are displayed on IDL plots.
  string base_unit;//Usually "cm" but sometimes "nm/s!u2!n"
  //Vector of vectors of doubles for the output.
  //1st index is the plot number- such as secular trend, period 1, period 2..
  //2nd index is the grid number.
  vector< vector<double> > outputs;

  //debug - records how many times each entry in output is changed.
  vector<long long> output_changed;

  //convert_results_from_ampphase_to_sincos() sets this value based on the
  //original type of results.
  //1 - Results were already in sincos.
  //2 - Results were in amp/phase and were converted to sincos.
  int conversion_to_sincos;

  //convert_results_from_sincos_to_ampphase() sets this value based on the
  //original type of results.
  //1 - Results were already in ampphase.
  //2 - Results were in sincos and were converted to ampphase.
  int conversion_to_ampphase;

  //Holds copy of options so I don't have to pass analysis_options around.
  analysis_options_s options;

  //Holds output for 2D XY plots.
  xy_s xy;

  //Holds output for equal-grid latlon maps.
  latlon_s latlon;

  //Some maps (e.g. type 4) have labels for each pt.
  vector< vector<string> > labels;

  //Each map now has an optional error bar. Any errors < 0 are ignored.
  vector<double> error_bars;
} results_s;

typedef struct {

  string outputfolder;//Folder where output is written.
  vector<string> output_files;//Filenames where outputs are written.
  string inputfolder;//Sometimes outputs are loaded, ex: load_output().

  int noplots;//Set to 1 if you don't want IDL to make plots.
  string subfolder;//Folder where other output is written.

  //write_univ_output creates a file with a list of filenames if this is the
  //final call of this fctn.  This file lists the output filenames so
  //IDL knows which files to open.  In order to insure that runs performed
  //on one computer are able to be plotted using IDL when running
  //on another computer, the filenames have to contain
  //ONLY the filename, not the path.  The path will be added by IDL later.
  vector <string> just_the_filenames;
} plot_options_s;

typedef struct {

  double lat,lon;
  vector<double> julian_times, sealevels;
  int type;//Controls how tides are loaded.
  int verbose;//Set to 1 to print predictions to screen.
} fes_data_s;

typedef struct {

  //If this is "0", have to calculate utc_seconds_in_year,jul_days_in_year
  //and choose dates to use in those vectors and populate "year" with them.
  //Also initialize dpm if this is 0.  Otherwise skip initialization steps.
  int initialized;

  //These 3 vectors are used to speed up multiple runs of conversion
  //functions by storing number of seconds in each year so that converting
  //multiple dates in the year 1936 to UTC seconds past J2000 only needs
  //to compute seconds in years 1937-1999 once.  Same with julian days.
  vector<long long> year;
  vector<long long> utc_seconds_in_year;
  vector<double> jul_days_in_year;

  //Stores indices of the years 1999 and 2000 in t_c.year for ease of use.
  long long pos1999,pos2000;

  //Holds days per month for every month.
  vector<int> dpm;

  //Holds calendar date/time.  Holds 6 values, here's what each index is:
  //0 - year
  //1 - month (1 = January)
  //2 - day (1 = first day in month)
  //3 - hour (0 = midnight,23=11pm)
  //4 - minute
  //5 - seconds
  vector<long long> calendar_time;

  //Holds UTC seconds since J2000:
  long long utc_seconds;

  //The FES Julian date is relative to jan-1-1950-00:00.00.
  //OR IS THAT OUTDATED?
  double fes_jul_days;

  double jd;//Standard Julian date.

  double dyr;//Decimal year (e.g. 2000.17).
} time_conversion_s;

typedef struct {

  //Used to record whether the span is negative, before absolute value is
  //taken for ywdhms routines.
  int negative;

  //Holds timespan in ywdhms format.  Holds 6 values:
  //0 - years
  //1 - weeks
  //2 - days
  //3 - hours
  //4 - minutes
  //5 - seconds
  vector<long long> ywdhms;

  long long seconds;//input
  double days;//input
  long long minute,hour,day,week,year;//holds number of seconds in each time unit
  vector<string> unit_titles;//Holds "year", etc for output string.
} timespan_conversion_s;

typedef struct {

  long long i,previous_i;//Current and previous index in the loop.
  long long num_loops;//How many loops total.
  long long interval;//The current index is modulo'd by interval for notifications
  long long time_remaining,time_loop,time_total;//Used for estimated time left.

  //Used to size notif. intervs, est time remaining
  time_t start_time,current_time,previous_time,end_time;

  //In seconds b/w notifications.
  long long min_interval, max_interval, first_interval;
  //notification_interval will not be allowed to go above this number
  //so that at LEAST a certain number of messages are seen.
  long long max_notif_interval;

  string prefix;

  //Used for converting timespans in seconds to human readable units.
  timespan_conversion_s timespan_conversion;

  //Next variable works the same as ps.perform_partial_saves, but it's
  //separate so that the ps version applies only to inversions,
  //while upward continuation can have ps.p_p_s == 1 and notif.p_s_e=0.
  //(Again, this is only because currently upward continuation doesn't
  //use partial saves. This may change in the future.)
  int partial_saves_enabled;

  //Set to 0 if not loading a partial save cache file, otherwise
  //set to ps.previous_input_index+1
  long long starting_index;
} notification_s;

typedef struct {

  int init;
  double sum,c,x;

  //Stores the size of vectors OR matrices.
  long long size;

  //Requires GSL libraries.
  #ifdef GSL_HERE
    //Used in ksum_vector()
    gsl_vector * sum_vector;
    gsl_vector * c_vector;
    gsl_vector * x_vector;
    gsl_vector * temp_vector1;

    //Used in ksum_matrix()
    gsl_matrix * sum_matrix;
    gsl_matrix * c_matrix;
    gsl_matrix * x_matrix;
    gsl_matrix * temp_matrix1;
  #endif
} kahan_s;

typedef struct {

  vector<double> a,b,c;//C = A X B
} cross_product_s;

typedef struct {

  double lat,lon,x,y,z,theta,costheta,sintheta;
  int init;
} transform_s;


//structure representing the astronomical angles.
typedef struct _astronomicAngle
{
    double  tt;
    double  hp;
    double  s;
    double  p1;
    double  p;
    double  iang;
    double  xi;
    double  nu;
    double  x1ra;
    double  r;
    double  nuprim;
    double  nusec;
} astronomicAngle;

//////////////////////////////////////////////////////////////////////////
//List of function declarations.
//////////////////////////////////////////////////////////////////////////

void create_options(string &options_file, int argc, char *argv[]);

void write_options(string &options_file,
                   input_s &input,
                   grid_options_s &region_options,
                   ancillary_s &ancillary,
                   grid_s &grid,
                   model_options_s &model_options,
                   analysis_options_s &analysis_options,
                   plot_options_s &plot_options,
                   insitu_s &insitu,
                   double parameter[],
                   long long parameter_type[],
                   long long num_parameters,
                   long long manual_list[],
                   long long num_manual_list,
                   double custom_lats[],
                   double custom_lons[],
                   long long num_custom_lats);

void load_options(string &input_file,
                  input_s &input,
                  grid_options_s &region_options,
                  ancillary_s &ancillary,
                  grid_s &grid,
                  model_options_s &model_options,
                  results_s &results,
                  plot_options_s &plot_options,
                  insitu_s &insitu,
                  int verbose);

void load_input(input_s &input, grid_s &grid, grid_s &region);

void clear_input(input_s &input);

void write_input(input_s &input);

input_s split_input(input_s &input, grid_s &region);

void load_ancillary(ancillary_s &ancillary, int verbose);

void define_grid(grid_s &grid, ancillary_s &ancillary, results_s &results);

void delete_grid_points_outside_boundaries(grid_s &grid);

void northern_latlons(grid_s &grid);

void southern_latlons(grid_s &grid);

void snap_grid_to_point(grid_s &grid);

void autosnap_grid_to_lat(grid_s &grid);

void snap_region_to_point(grid_s &region);

long long index_of_closest_match(vector<double> &dbls, double &temp);

void calculate_grid_xyz(grid_s &grid);

void xyz2llh(double x, double y, double z,
             double &lat, double &lon, double &h);

void llh2xyz(double &x, double &y, double &z,
             double lat, double lon, double h);

void calculate_grid_topography(grid_s &grid, ancillary_s &ancillary);

void calculate_grid_topography2(grid_s &grid, ancillary_s &ancillary);

void rotate_grid_about_x(grid_s &grid);

void rotate_grid_about_z(grid_s &grid);

void autosize_grid_cap(grid_s &grid);

void autosize_region_cap(grid_s &region);

void grid_cm2kg(grid_s &grid);

void write_grid_to_cache(grid_s &grid);

void load_grid_from_cache(grid_s &grid);

void define_support_grid(grid_s &grid);

void insitu_initialize(insitu_s &insitu,
                       grid_s &grid,
                       model_options_s &model_options,
                       ancillary_s &ancillary);

void load_insitu(insitu_s &insitu, ancillary_s &ancillary);

void load_amps_n_phases(insitu_s &insitu, long long &current_station_index);

void load_proshutinsky(insitu_s &insitu, long long &current_station_index);

void adjust_tg_phases(insitu_s &insitu, long long index);

void generate_tg_timeseries(insitu_s &insitu, long long index);

grid_s define_region(grid_options_s &region_options,
                     grid_s &grid,
                     int verbose);

void search_for_grid_pts_in_region(grid_s &region,
                                   grid_s &grid,
                                   int verbose);

void calculate_region_as_grid(grid_s &region);

void load_model(model_s &model, grid_s &grid);

void define_model_titles_and_periods(model_s &model);

void read_rawfes(latlon_s &rawfes);

void write_rawfes(latlon_s &rawfes);

void add_residuals_to_fes(results_s &results, grid_s &grid, model_s &model);

void match_model_and_results_parameters(results_s &results, model_s &model);

void build_analysis_region_list(grid_s &grid, grid_s &region);

void perform_analysis(input_s &input,
                      grid_s &grid,
                      ancillary_s &ancillary,
                      model_s &model,
                      results_s &results,
                      plot_options_s &plot_options,
                      grid_s &region,
                      insitu_s &insitu);

void perform_timeseries_analysis(input_s &input, results_s &results);

void sequential_regional_analysis(
                     input_s &input,
                     grid_s &grid,
                     ancillary_s &ancillary,
                     model_s &model,
                     results_s &results,
                     results_s &regional_results,
                     results_s &overlap_results,
                     plot_options_s &plot_options,
                     plot_options_s &temp_options,
                     plot_options_s &overlap_plot_options,
                     grid_s &region,
                     grid_s &empty_region,
                     insitu_s &insitu);

void parallel_regional_analysis(
                     input_s &input,
                     grid_s &grid,
                     ancillary_s &ancillary,
                     model_s &model,
                     results_s &results,
                     results_s &regional_results,
                     results_s &overlap_results,
                     plot_options_s &plot_options,
                     plot_options_s &temp_options,
                     plot_options_s &overlap_plot_options,
                     grid_s &region,
                     grid_s &empty_region,
                     insitu_s &insitu);

//Requires GSL libraries.
#ifdef GSL_HERE
  void copy_x_to_outputs(results_s &results, grid_s &grid, gsl_vector * x);
  
  void copy_x_to_ampphase(gsl_vector * x,
                          results_s &ar_series1,
                          insitu_s &insitu,
                          long long i);
#endif

void ksum(kahan_s &kahan);

void ksum_vector(kahan_s &kahan);

void ksum_matrix(kahan_s &kahan);

h_matrix_s create_h_matrix(input_s &input, grid_s &grid, results_s &results);

least_squares_s create_least_squares_in_parallel(input_s &input,
                                                 grid_s &grid,
                                                 results_s &results);

void least_squares_load_from_cache(least_squares_s &least_squares,
                                   results_s &results,
                                   grid_s &grid);

void least_squares_save_to_cache(least_squares_s &least_squares,
                                 results_s &results,
                                 input_s &input);

void mirror_across_diagonal(least_squares_s &least_squares,
                            results_s &results);

void define_zone_of_influence(grid_s &grid, results_s &results);

void scale_b_vector(least_squares_s &least_squares,
                    results_s &results,
                    double scale_factor);

void allocate_least_squares_matrices(
                          least_squares_s &least_squares,
                          grid_s &grid,
                          results_s &results);

void create_least_squares(least_squares_s &least_squares,
                          input_s &input,
                          grid_s &grid,
                          results_s &results);

void notify(notification_s &notification, long long i,
            least_squares_s &least_squares,
            input_s &input,
            results_s &results);

void least_squares_cache_initialize(grid_s &grid,
                                    grid_s &region,
                                    input_s &input,
                                    results_s &results);

void read_partial_save_details_file(results_s &results);

void write_partial_save_details_file(results_s &results);

void check_ls_subfolder_completeness(results_s &results, grid_s &region);

void perform_partial_save(least_squares_s &least_squares,
                          input_s &input,
                          results_s &results,
                          long long &i);

least_squares_s regularize(grid_s &grid, results_s &results);

void regularize_in_place(grid_s &grid,
                         results_s &results,
                         least_squares_s &least_squares);

void define_regional_grid(grid_s &grid,
                          grid_s &region,
                          grid_s &regional_grid,
                          long long &i);

input_s load_regional_input(long long regionindex,
                            input_s &maininput,
                            grid_s &region);

void copy_regional_results(results_s &regional_results,
                           results_s &results,
                           grid_s &region,
                           grid_s &grid,
                           grid_s &regional_grid,
                           long long &i);

void display_output(results_s &results, grid_s &grid, long long number);

void write_output(string base,
                  results_s &results,
                  grid_s &grid,
                  plot_options_s &plot_options,
                  int final);

void init_filenames(string base,
                    long long numfiles,
                    plot_options_s &plot_options);

void write_map_header(results_s &results, long long j, FILE *new_fp);

void create_output_filenames(plot_options_s &plot_options,
                             results_s &results);

void find_and_load_parallel_region_output(
              grid_s &region,
              plot_options_s &plot_options,
              results_s &results,
              grid_s &grid);

void read_output(results_s &results,
                 string filename,
                 grid_s &grid);

void insitu_comparison_verification(insitu_s &insitu,
                                    results_s &results,
                                    grid_s &grid,
                                    model_options_s &model_options,
                                    plot_options_s &plot_options);

void upward_continuation_verification(grid_s &grid,
                                      input_s &input,
                                      input_s &unused_input,
                                      results_s &results,
                                      grid_s &region,
                                      plot_options_s &plot_options);

void copy_plot_line(results_s &temp_results,
                    results_s &noise,
                    long long plotnum,
                    long long linenum,
                    long long newplot);

void extract_signal(input_s &input, long long verif, long long type);

void create_synth_timeseries(insitu_s &insitu);

void generate_model_timeseries(insitu_s &insitu, model_options_s &model_options);

void add_residuals_to_fes(results_s &results, ancillary_s &ancillary, grid_s &grid);

void generate_input_timeseries(insitu_s &insitu,
                               results_s &results,
                               grid_s &grid,
                               int subtract_mean);

void compare_station_timeseries(insitu_s &insitu);

results_s analyze_timeseries(insitu_s &insitu);

void make_timeseries_plots(insitu_s &insitu);

double spherical_distance(double lat1, double lon1, double lat2, double lon2);

void spherical_latlon2xyz(transform_s &transform);

void spherical_xyz2latlon(transform_s &transform);

void rotate_xyz_about_x(transform_s &transform);

string create_grid_cache_filename(grid_options_s &grid_options);

string create_grid_options_substring(grid_options_s &grid_options);

string create_support_grid_cache_filename(grid_options_s &grid_options);

string create_model_cache_filename(grid_s &grid, model_s &model);

string create_region_cache_filename(grid_options_s &region_options,
                                    grid_options_s &grid_options);

string create_region_options_substring(grid_options_s &region_options);

string create_least_squares_cache_foldername(
                      grid_options_s &region_options,
                      grid_options_s &grid_options,
                      results_s &results);

string create_input_cache_filename(input_options_s &input_options,
                                   grid_options_s &grid_options);

string create_tg_cache_filename(insitu_s &insitu);

string create_least_squares_cache_filename(results_s &results);

string create_arbitrary_cache_filename(input_s &input);

void swap_double(double *i, double *j);

void swap_long(long long *i, long long *j);

int is_pt_in_bin(grid_s &grid, double current_lat, double current_lon,
                 long long i, int &newbin);

int is_pt_in_grid(grid_s &grid, double current_lat, double current_lon);

void wrap_longitudes(grid_s &grid);

long long grid_index(grid_s &grid, double current_lat,
                     double current_lon, int verbose);

long long grid_index_via_distance(grid_s &grid, double current_lat,
                                  double current_lon, int verbose);

double nearest_topo_pt_elev(double lat, double lon,
                            ancillary_s &ancillary);

grid_s local_topography(double current_lat, double current_lon,
                        ancillary_s &ancillary);

grid_s local_topography2(grid_s &grid,
                         long long i,
                         ancillary_s &ancillary);

grid_s load_fes_points_near_insitu_stations(insitu_s &insitu,
                                            ancillary_s &ancillary,
                                            model_s &model);

void search_input_cache_for_binindices(input_s &input,
                                       string &filename,
                                       int &matching_cache_found);

void compare_input_files(input_s &input);

void compare_input_file_to_verification_file(input_s &input,
                                             grid_s &grid,
                                             results_s &results);

//I did not originally write readdoubles- it was copied from the
//FES2004 source code as "readValues".  I simply modified it slightly
//to read doubles.  At the moment I don't really understand
//why its return value is a "static int" or why its inputs are
//"const int" and "double*" (though double* is probably a pointer
//which is how value is able to be used to return... um... values.
int readdoubles(char* line, const int sizeOf,double* values);

void astronomics(const double itj, astronomicAngle* angle);

int fes_access(fes_data_s &fes_data, insitu_s &insitu);

double sqr(const double x);

void create_output_vectors(results_s &results, input_s &input);

void greenwich_phase_at_epoch(doodson_s &doodson);

void cp_doodson(int doodson_numbers[], period_info_s &period_info);

void is_this_period_recognized(double period, period_info_s &period_info);

void convert_title_to_parameter(string title, period_info_s &period_info);

void init_time_conversion(time_conversion_s &time_conversion);

void cal2utc(time_conversion_s &time_conversion);

void dyr2utc(time_conversion_s &time_conversion);

void utc2dyr(time_conversion_s &time_conversion);

void utc2cal(time_conversion_s &time_conversion);

void cal2fes_jul(time_conversion_s &time_conversion);

void utc2fes_jul(time_conversion_s &time_conversion);

void jd2utc(time_conversion_s &time_conversion);

string sec2human(timespan_conversion_s &timespan_conversion);

string sec2cal(long long i);

string seconds_to_string(long long i);

string day2human(timespan_conversion_s &timespan_conversion);

string current_datetime();

void cross(cross_product_s &cross_product);

void load_output(results_s &results,
                 plot_options_s &plot_options,
                 grid_s &grid);

void convert_model_to_results(grid_s &grid, results_s &results, model_s &model);

void convert_rawfes_to_latlon(results_s &results);

void convert_latlon_to_rawfes(results_s &results);

void convert_results_from_ampphase_to_sincos(results_s &results);

void convert_results_from_sincos_to_ampphase(results_s &results);

void convert_model_from_ampphase_to_sincos(model_s &model, grid_s &grid);

void convert_model_from_sincos_to_ampphase(model_s &model, grid_s &grid);

void adjust_phases(input_s &input, results_s &results);

double restrict_degrees(double phase);

double restrict_degrees360(double phase);

double restrict_radians(double phase);

void create_folder(string &folder);

void remove_folder(string &folder);

string gaiafolder();

void create_output_folder(int argc, char *argv[], plot_options_s &plot_options);

//These functions are used in combine_orbit_and_accels.cpp

void append_input_and_orbit_data(input_s &input, long long &mainindex);

input_s load_another_input_segment(input_options_s &input_options,
                                   long long &mainindex);

//These functions are used in sim_accels.cpp

void perform_upward_continuation(input_s &input, grid_s &grid, results_s &results);

void perform_upward_continuation2(input_s &input,
                                  grid_s &grid,
                                  results_s &results,
                                  results_s &ntresults,
                                  long long accelmapindex);

void save_regional_input(string regioninputfilename, input_s &input);

void save_regional_verification_input(string regioninputfilename,
                                      input_s &input);

void combine_regional_inputs(input_s &input, grid_s &region);

void combine_regional_verification_inputs(input_s &input, grid_s &region);

void copy_results_from_global_to_regional(
                         results_s &regional_results,
                         results_s &results,
                         grid_s &region,
                         grid_s &grid,
                         grid_s &regional_grid,
                         long long &i);

void copy_output_to_latlon(results_s &results, grid_s &grid);

void copy_results_to_latlon(results_s &results,
                            grid_s &grid,
                            plot_options_s &plot_options);

void write_supporting_grid(grid_s &grid);

void estimate_land_tide(results_s &results, grid_s &grid);

void relative_reduction(results_s &results);

void subtract_fes(latlon_s &rawfes, latlon_s &rawfes2, int mode);

void decimate_fes(latlon_s &rawfes,
                  latlon_s &rawfes2,
                  int decimate);

void parse_input_line(char s[max_length], long long &onetime, double &onelat, double &onelon, double &oneaccel, int &oneascending, double &onexA, double &onexB, double &oneyA, double &oneyB, double &onezA, double &onezB, int no_orbit_data);

void match_fes_to_grid(latlon_s &rawfes,grid_s &grid);

void sum_fes_areas_in_each_mascon(grid_s &grid);

void scale_plot(results_s &results,int target,int source);//Scale line at first index by lines in second index (-1 means all others)

void combine_grids(grid_s &grid1,grid_s &grid2, int choice);

void load_arbitrary_time_series(input_s &input, results_s &results);

void read_arbitrary_file(xy_s &arbitrary, long long type);

void arbitrary_interpolate(input_s &input, grid_s &region);

void interpolate_preloaded_input_and_arbitrary(input_s &input);

void closest_arbitrary_times(input_s &input, long long i, long long j);

int is_this_arbitrary_point_closer(input_s &input, long long i,
                                   long long j, long long k);

void save_interpolated_arbitrary_to_cache(input_s &input);

//This is the include guard.
#endif
